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ABSTRACT 


a 

■^Analytical  response  predictions  (pre-shot)  are  made  for  three  structures 
that  are  to  be  evaluated  during  the  Mighty  Epic  tests.  Also  included 
in  the  report  is  a prediction  of  the  response  for  the  free  field  and  for 
the  free  field  containing  a cavity,  as  well  as  a study  of  the  effects  that 
various  idealization  parameters  have  on  the  predicted  results. > 

OBJECTIVE  AND  SCOPE 

The  objective  of  this  study  is  to  analytically  predict  (pre-shot)  the 
response  of  three  structures  that  are  in  place  at  the  Mighty  Epic  test 
site.  These  structures,- shown  in  Figures  1 and  consist  of: 

;^ll  Homogeneous  sphere,  ^ 

Composite  built-up  liner^ 

3)  Composite  integral  liner. 

These  three  structures  are  situated  in  a relatively  weak  rock  (tuff)  at 
the  0.5-Kbar  range.  Included  in  the  study  is  a prediction  of  the  response 
for  the  free  field  (i.e.,  tuff)  and  for  the  tuff  containing  a six-foot 
diameter  cylindrical  cavity.  Also  reported  are  the  effects  that  the  various 
parameters  associated  with  the  predictions  have  on  the  results. 

The  scope  of  the  analyses  is  predicated  on  usage  of  NONSAPl  computer 
code  in  its  "as  is"  state.  NONSAP  is  a finite  element  computer  program 
designed  primarily  for  the  nonlinear  analysis  of  three-dimensionalx  structures 
subjected  to  either  static  or  dynamic  forces  applied  either  separately 
or  in  combination.  The  program  offers  a variety  of  material  models,  two 
implicit  time  integration  schemes  (i.e.,  Wilson's  6 or  Newmark's  B method), 
two  finite  deformation  formulations  (based  on  either  an  Eulerian  or  Lagrangian 
reference  system),  and  several  element  types  (including  the  truss,  beam, 

4-8  node  isoparametric  quadrilateral,  8-21  node  isoparametric  hexahediron, 
and  several  types  of  joint  and  boundary  elements).  An  "in-core"  version 
of  NONSAP,  operated  on  a CDC  7600,  is  used  to  develop  the  results  presented 
in  this  report. 

This  report  was  prepared  for  use  as  an  internal  working  document. 
Therefore,  its  distribution  is  limited,  and  disclosure  of  all  or  part  of 
its  contents  outside  the  Government  must  have  prior  approval  of  the  Civil 
Engineering  Laboratory  or  the  sponsor  of  the  work  reported. 
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PROBLEM  DESCRIPTION 
NONSAP  Parameters 

NONSAP  is  particularly  designed  to  solve  nonlinear  problems.  The 
basic  algorithm  uses  an  incremental  tangent  approximation  for  the  nonlinear 
stiffness  coupled  with  equilibrium  iterations  at  the  end  of  a step  to  correct 
for  the  errors  which  arise  as  a consequence  of  the  tangent  approximation. 

For  the  problems  discussed  in  this  report,  the  stiffness  matrix  is  inverted 
at  every  step  while  equilibriimi  iteration  is  used  only  for  problems  where 
large  deformations  are  included.  The  large  deformation  formulation  is 
based  on  a Lagrangian  reference  system  which  is  felt  to  be  most  compatible 
with  the  material  data  provided.  However,  the  lack  of  large  deformation- 
oriented  material  data  debilitates  the  employment  of  this  formulation. 

The  time  integration  is  accomplished  using  Wilson's  6 method. 

For  the  problems  in  this  report,  16  steps  are  used  in  application  of 
nonlinear  static  runs  while  44  steps  with  a 0.4-msec  time  step  are  used 
for  dynamic  runs  (providing  17.6  msec  of  response).  Typical  CDC  7600  central 
processing  times  and  costs  are: 

o For  the  nonlinear  static  cases  with  equilibrium  iteration, 

210  seconds  amd  $80. 

o For  the  linear  static  cases,  7 seconds  and  $10. 

o For  the  nonlinear  dynamic  cases  with  equilibrium  iteration, 

470  seconds  and  $160. 

o For  the  linear  dynamic  cases,  45  seconds  and  $15. 

The  high  cost  of  equilibrium  iteration  is  caused  by  its  inefficient 
iii5>lementation  on  the  CDC  7600.  Reprogramming  this  algorithm  in  accordance 
with  the  7600  architecture  would  mitigate  this  problem. 

The  elements  employed  in  all  solutions  are  4-8  node,  isoparametric 
quadrilaterals.  Integration  orders  from  2-4  may  be  specified;  in  general, 
orders  of  3 are  used  for  the  computations  in  this  report.  Integration 
orders  of  3 result  in  nine  integration  points  being  established  within 
the  element  at  the  Gauss  points.  Separate  stress  histories  and  material 
properties  are  computed  for  each  of  these  points. 

Some  additional  relevant  features  of  NONSAP  are: 

1.  Its  restart  capability,  whereby  any  solution  may  be  stopped  and 
restarted  with  different  parameters.  For  example,  the  geostatic  stress 
at  the  Mighty  Epic  site  may  be  computed  as  a nonlinear  statics  problem; 
these  results  would  be  automatically  stored  and  be  available  as  initial 
boundary  conditions  for  the  dynamic  NONSAP  simulation  of  the  Mighty  Epic 
structures. 

2.  Birth  and  death  of  elements,  whereby  elements  are  assigned  times 
of  appearance  and/or  disappearance.  This  makes  it  possible  to  solve  the 
Mighty  Epic  geostatic  problem  as  an  xinlined  cavity  where  prior  to  the 
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dynamic  run,  the  structure  and  grout  elements  would  be  "bom"  into  the 
structural  system. 

General  Mesh  Description 

Five  different  structural  problems  are  presented  in  this  study: 

1.  Free  field 

2.  Free  field  with  cavity 

3.  Homogeneous  sphere 

4.  Composite  built-up  liner 

5.  Composite  integral  liner 

Each  structure  was  analyzed  in  four  basic  configurations:  (1)  linear 
static,  (2)  nonlinear  static,  (3)  linear  dynamic,  and  (4)  nonlinear  dynamic. 
Within  each  of  the  four  basic  configurations,  several  parameters  were 
varied  (e.g.,  mesh  size,  time  step,  degree  and  type  of  nonlinearity)  with 
the  intent  of  determining  their  effect  and  developing  confidence  in  the 
final  prediction  for  the  nonlinear  dynamic  response. 

The  general  character  of  the  finite  element  mesh  used  for  the  five 
structural  problems  is  shown  in  Figure  3.  The  column  idealization  is  used 
for  the  free  field  calculations,  while  the  cavity  idealization  is  used 
for  the  other  four  structures.  The  center  of  the  cavity  is  set  252  inches 
(approximately  6 meters)  below  the  surface  of  the  mesh.  The  side  boundaries 
of  both  meshes  are  constrained  horizontally  while  the  mesh  bottom  is  fixed 
far  enough  from  the  region  of  interest  to  prevent  any  interference  caused 
by  reflecting  waves.  The  pressure  loads  are  applied  uniformly  across  the 
top  of  the  structure  either  statically  or  as  functions  of  time  (dynamic). 
Typically,  a thousand  degrees-of-freedom  are  used  to  model  the  cavity  ideal- 
izations while  fewer  are  needed  for  the  free  field  modeling. 

Materials 

These  structures  are  composed  of  five  different  types  of  materials: 

Tuff 

Steel 

Fiber  reinforced  concrete 

Reinforced  concrete 

Cellular  concrete 

The  material  constants  are  listed  in  Table  1;  their  derivation  is 
given  in  Appendix  A.  Two  sets  of  properties  are  given  for  tuff:  Tuff 
No.  2 is  based  on  in-situ  testing  and  Tuff  No.  1 represents  an  earlier 
assessment  of  material  behavior  generated  prior  to  the  in-situ  testing. 

Geostatic  Effects 

The  gravity  effects  are  complicated  by  the  employment  of  tunneling 
and  grouting  for  structure  emplacement.  NONSAP  could  compute  the  geostatic 
stresses  in  tuff  surrounding  the  tunnel.  However,  since  the  grout 
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surrounding  the  structures  is  in  an  unstressed  condition  (at  least  from 
gravity),  and  its  dimensions  uncertain  (possibly  2-4  feet  between  structure 
and  tunnel,  and  in  the  case  of  the  sphere  essentially  unlimited  along  the 
tunnel  axis),  the  inclusion  of  geostatic  stresses  seems  unwarranted  and 
somewhat  specious.  Also,  since  the  estimated  geostatic  stress  in  the  free 
field  at  the  emplacement  depth  is  only  70  bars^,  it  does  not  appear  the 
omission  of  gravitational  forces  will  produce  a serious  degradation  in 
; predicted  response. 

j Forcing  Function 

r 

I The  forcing  functions  considered  in-this  study  are  derived  from  the 

; curves  (supplied  by  Weidlinger  and  Associates^)  shown  in  Figures  4 and 

5.  They  consist  of  uniaxial  velocity  and  stress  at  various  locations  with 
respect  to  the  structure  (where,  for  example,  -6m  denotes  a point  that 
j is  six  meters  upstream  from  the  structure's  center).  These  curves  are 

derived  analytically  based  on  a CAP  model'*  characterization  of  the  tuff. 

( From  these  curves  and  through  a study  of  the  free  field  calculations  (next 

' section),  the  forcing  function  for  the  dynamic  configurations  was  chosen. 

TVo  forms  of  this  function  are  shown  in  Figure  6:  one  is  a detailed  descrip- 
tion of  the  curve  which  was  applied  to  all  five  structural  problems  while 
the  other  is  a simplified  version  developed  to  aid  in  studying  the  effects 
• of  various  parameters  on  the  free  field  response.  For  the  static  analyses, 

i a peak  pressure  of  7250  psi  was  used. 

FREE  FIELD  VALIDATION  AND  DETERMINATION  OF  VARIOUS  MESH  PARAMETERS 

Free  Field  Validation 

This  section  presents  the  studies  conducted  to  validate  the  accuracy 
of  the  free  field  responses  computed  by  NONSAP.  The  basic  premise  is  that 
the  free  field  response  which  occurs  in  the  tuff/structure  meshes  (presented 
in  next  section)  is  similar  to  what  occurs  in  a columnar  mesh  (Figure  7) 
made  of  tuff  only.  The  element  and  step  size  selected  for  the  columnar 
mesh  are  based  on  the  results  presented  in  ^pendix  B.  The  elements  are 
20  inches  in  height,  4-node  quadrilaterals,  constrained  laterally,  and 
subjected  at  the  mesh  surface  to  the  -6m  stress  wave  shown  in  Figure  4. 

A time  step  of  0.4  msec  is  used  to  compute  the  results  shown  in  Figures 
8 and  9.  Figure  8 shows  stress  versus  time  at  element  24  for  the  linear 
elastic  and  the  Drucker-Prager  approximations  of  Tuff  No.  1.  Figure  9 
shows  velocity  versus  time  for  node  49.  The  +6  meter  CAP  model  results, 
which  correspond  in  depth  to  element  24  and  node  49,  are  also  shown  in 
the  figures,  as  well  as  the  stress  wave  applied  to  the  mesh. 

The  results  of  these  free  field  runs  indicate: 

1.  Comparisons  of  the  vertical  stress/velocity  is  not  an  appropriate 
method  for  checking  the  validity  of  the  tuff  material  characterization. 

More  germane  would  be  the  comparison  of  horizontal  stresses  or  material 
stiffness  predictions  at  the  springline/crown  of  a lined  cavity.  These 
stresses  are  much  more  sensitive  to  variations  between  the  various  material 
models. 
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2.  Velocity  output  does  not  compare  as  well  as  stress  output. 
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3.  No  appreciable  decrease  in  the  peak  stress  with  travel  distance 
occurs  for  either  the  linear  or  Drucker-Prager  models.  Moreover,  these 
two  models  yield  similar  results  for  the  columnar  idealization. 

Compatibility  of  4-Node  and  Higher  Order  Quadrilateral  Elements 

In  certain  regions  of  a mesh  it  is  advantageous  to  join  two  4 node 
quadrilateral  elements  to  the  side  of  an  element  of  higher  order  (i.e., 
one  with  5 to  8 nodes).  Figure  10a  shows  such  a case.  Because  the  midside 
node  of  a higher  order  element  is  4 times  stiffer  than  its  comer  node, 
some  anomaly  is  expected  as  a stress  wave  passes  through  the  transition 
region.  To  investigate  this  phenomenon,  the  top  surfaces  of  the  three 
free  field  meshes  shown  in  Figure  10  were  subjected  to  a ramp  loading  of 
infinite  duration.  The  results  of  this  study  are  given  in  Appendix  C. 

These  results  indicate  that  while  a local  anomaly  in  the  displacement/stress 
field  transpires,  the  effect  at  a point  some  distance  from  this  region 
is  negligible.  Therefore,  these  transition  zones  will  not  be  located  near 
regions  of  interest. 

Technique  Employed  for  Removal  of  Reflection  Waves 

A significant  aspect  of  predicting  the  transient  response  is  insuring 
that  the  fictitious  boundaries  within  the  free  field  imposed  by  the  finite 
element  idealization  do  not  become  sources  for  spurious  reflection  waves. 
Because  the  forcing  functions  are  uniaxial,  the  location  and  type  of  side 
boundaries  employed  appears  to  have  little  effect  on  the  computations  (see 
results  for  Mesh  3.0  and  3.1  of  next  section).  Reflections  caused  by  the 
bottom  boundary  can  most  simply  be  dealt  with  by  establishing  its  location 
far  enough  from  the  region  of  interest.  The  depth  of  the  region  of  interest 
is  denoted  as  Dj  while  Db  is  the  additional  depth  of  mesh  necessary  to 
avoid  reflection  waves.  Figure  3.  Dg  is  computed  from  the  wave  speed  in 
tuff  (C  = 57.5  in/msec,  from  Appendix  B)  and  the  amount  of  time  needed 
to  capture  the  response  (estimated  to  be  no  more  than  35  msec) . Therefore, 
the  depth  of  mesh  needed  beyond  the  region  of  interest  is: 

Dg  = (57.5)  (35) /2  * 1000  inches 

Selecting  the  depth  of  the  region  of  interest  (Dj)  to  be  500  inches,  yields 
the  total  depth  of  mesh: 

D = Dg  + Dj  = 1500  inches 

Since  the  element  size  needed  in  this  region  is  40  inches  (see  Appendix  B) 
the  total  number  of  degrees  of  freedom  needed  to  mcdel  this  space  is: 

DOF  (for  4-node  quads)  = 1000/40  (2)  = 50 

DOF  (for  8-node  quads)  = 1000/40  (5)  = 125 
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where  one  element  is  presumed  to  cover  the  whole  space  laterally.  Whichever 
element  type  is  employed,  the  niunber  of  degrees-of-freedom  used  to  model 
this  space  is  inconsequential  compared  to  the  total  degrees-of-freedom 
for  the  whole  model  (i.e.,  more  than  1000).  Using  more  complicated  schemes, 
such  as  nodal  dampers,  does  not  seem  warranted  in  these  circumstances. 

Computation  of  Residual  Response 

A possible  scheme  for  computing  residual  response  consists  of  a deep 
mesh  of  successively  larger  element  sizes  used  in  conjunction  with  the 
NONSAP  restart  option,  whereby  a small  time  step  is  used  until  the  peak 
response  of  the  structure  occurs . Then  gradually  larger  time  steps  are 
employed  until  the  structure  is  in  static  equilibrium.  An  example  of  this 
technique  is  shown  in  Appendix  D.  In  general,  for  these  types  of  problems, 
this  technique  appears  to  provide  a way  to  reduce  computer  costs  for  large 
runs  and  to  be  a feasible  means  of  computing  response  over  a long  period 
of  time. 

Use  of  Velocity  as  Forcing  Function 

NONSAP  does  not  allow  direct  input  of  velocity  as  a forcing  function. 

To  use  velocity  data,  it  must  be  integrated  to  provide  a displacement-time 
history.  Then  using  the  "stiff  spring  trick"  (described  in  Appendix  E), 
these  displacements  are  used  to  drive  the  problem.  Since  both  pressure 
and  vel  data  were  provided,  there  was  no  need  to  perform  these  operations 

for  a:  ne  analyses  reported.  However,  for  comparative  purposes,  a 

free  culation  based  on  these  techniques  is  included  in  ^pendix  E. 

ANA  Jt-'  A CYLINDRICAL  CAVITY  WITHIN  THE  FREE  FIELD 

Four  meshes  were  made  to  model  a cylindrical  cavity  (72  inches  in  dia- 
meter) within  a homogeneous  tuff.  Such  parameters  as  mesh  boundaries; 
element  sizes,  types,  and  shapes;  and  linear  versus  nonlinear  materials 
and  deformation  were  studies.  These  studies  were  carried  out  for  both 
an  unlined  cavity  and  one  lined  with  a thin  steel  shell  (0.78  inches  thick). 
The  steel- lined  cavity  results  are  included  to  make  the  parametric  studies 
relevant  to  the  Mighty  Epic  structures. 

For  this  analysis,  seven  static  and  seven  dynamic  problem  configurations 
were  deemed  to  be  of  interest.  These  are: 

1.  An  elastic  tuff  with  an  unlined  cavity  (denoted  SCI  and  DCl). 

2.  A Drucker-Prager  tuff  with  an  unlined  cavity  (denoted  SC2  and  DC2). 

3.  An  elastic  tuff  with  an  unlined  cavity  using  the  large  deformation 
formulation  and  equilibrium  iteration  (denoted  SC3  and  DC3). 

4.  A Drucker-Prager  tuff  with  an  unlined  cavity  using  the  large  deform- 
ation formulation  and  equilibrium  iteration  (denoted  SC4  and  DC4) . 

5.  An  elastic  tuff  with  an  elastic  steel-lined  cavity  (denoted  SC5 
and  DCS) . 
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6.  A Drucker-Prager  tuff  with  an  elastic  steel-lined  cavity  (denoted 
SC6  and  DC6) . 

7.  A Drucker-Prager  tuff  with  an  elastic  steel-lined  cavity  using 
the  large  deformation  formulation  and  equilibrium  iteration  (denoted  SC7 
and  DC7) . 

The  following  four  meshes  were  used: 

Mesh  3.0:  Conventional  4-node  quadrilateral  mesh  where  thj  right 
boundary  is  set  at  Y = 288  inches.  The  thin  liner  consists  of  thirty- two 
8-node  quadrilaterals  connected  by  5-node  quadrilaterals  to  the  rest  of 
the  mesh.  Figure  11a. 

Mesh  3.1:  Identical  to  Mesh  3.0  except  that  additional  elements  have 
been  added  to  the  right  side  moving  the  right  boundary  to  Y = 435  inches. 
Figure  11b. 

Mesh  3.2:  Similar  to  Mesh  3.0  except  that  8-node  quadrilaterals  have 
(approximately)  replaced  groups  of  four  4-node  quadrilaterals.  The  thin 
liner  consists  of  only  sixteen  8-node  quadrilaterals.  Figure  11c. 

Mesh  3.3:  Provides  a rectangular  mesh  for  contrast  with  the  previous 
radial  meshes.  The  thin  liner  consists  of  eighteen  8-node  quadrilaterals 
which  are  connected  by  5-node  quadrilaterals  to  the  rest  of  the  4-node 
quadrilateral  mesh. 

Only  the  upper  region  of  each  mesh  is  shown  in  Figure  11.  The  lower  re- 
gions are  similar  in  depth  to  Figure  Sb. 

These  four  meshes  are  presumed  to  provide  sufficient  data  for  evaluation 
of  right  boundary  selection;  element  sizing,  layout  and  type;  and  the  effect 
of  mesh  selection  on  response  predictions.  Both  static  and  dynamic  responses 
were  obtained  at  the  mesh  locations  shown  in  Figure  12.  The  static  responses 
were  included  in  Table  2-5,  while  the  peak  dynamic  responses  were  included 
in  Tables  6-9.  A typical  dynamic  response  for  the  four  meshes  is  shown 
in  Figure  13.  The  seven  problem  configurations  were  not  run  for  every 
mesh,  owing  to  expected  similarities  of  results. 

The  response  of  each  mesh  for  a particular  problem  configuration  is 
essentially  the  same.  Comparison  of  Mesh  3.3  to  the  other  meshes  indicates 
that  the  easily  generated  radial  meshes  can  satisfactorily  be  employed 
in  the  structural  meshes.  Comparison  of  Meshes  3.0  and  3.1  indicates  that 
a right  boundary  of  288  inches  is  satisfactory.  Comparison  of  Meshes  3.0 
and  3.2  indicates  that  sixteen  8-node  quadrilaterals  are  adequate  to  handle 
the  bending  present  in  a thin  liner,  and  8-node  quadrilaterals  can  replace 
groups  of  four  4-node  quadrilaterals  in  the  medium.  Comparison  of  Runs 
SC6  and  SC7  indicates  that  when  the  steel  liners  are  employed,  deflections 
are  so  small  that  large  deformation  theory  is  not  necessary. 
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RESULTS  OF  ANALYSES  OF  MIGHTY  EPIC  STRUCTURES 


Based  on  the  results  of  the  previous  section  and  various  finite  element 
analyses  of  the  structures  alone,  the  following  three  structural  meshes 
were  generated; 

Mesh  4.0:  Homogeneous  Sphere  - Conventional  4-node  quadrilateral  mesh 
where  the  right  boundary  is  set  at  Y = 288  inches.  The  concrete  sphere 
consists  of  six  layers  of  4-node  quadrilaterals. 

Mesh  5.0;  Composite  Built-Up  Liner  - Conventional  8-node  quadrilateral 
mesh.  The  steel  liner  and  cellular  concrete  consists  of  1 and  2 element 
layers,  respectively. 

Mesh  6.0;  Composite  Integral  Liner  - Identical  in  topology  to  Mesh  5.0 

The  upper  region  of  each  mesh  is  shown  in  Figures  14,  15  and  16.  For  each 
mesh,  four  different  runs  were  made: 

1.  Linear  static. 

2.  Nonlinear  static,  where  the  nonlinear  material  properties  of  Table 
1 were  employed  and  small  deformation  theory  was  assumed  valid. 

3.  Linear  dynamic,  where  At  = 0.4  msec  and  44  steps  were  employed. 

4.  Nonlinear  dynamic,  where  the  material  characterizations  of  Run  2 
and  the  time  parameters  of  Run  3 were  employed. 

Homogeneous  Sphere  (Figure  1) 

Since  the  designers  of  the  homogeneous  sphere  based  their  design  on 
the  stress  distribution  within  a thick-walled  sphere  subjected  to  a uniform 
external  pressure,  the  elastic  thick-wall  theory  in  its  general  form  is 
given  below: 


where : 

= wall  stress  in  tangential  direction  (psi) 

Sy  = wall  stress  in  radial  direction  (psi) 

Po  = external  pressure  (psi) 
ro  = external  radius  (in) 
rj  = interior  radius  (in) 

r = radius  to  location  in  wall  under  consideration  (in) 
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Using  these  equations,  an  external  pressure  of  7250  psi  acting  on  the  sphere 
shown  in  Figure  1 would  produce  tangential  and  radial  stresses  of  -15,450 
psi  and  0 psi  at  the  inner  surface  and  -11,820  psi  and  -7250  psi  at  the 
outer  surface.  The  diameter  change  of  the  sphere  with  the  elastic  proper- 
ties shown  in  Table  1 would  equal  -0.17  in. 

The  static  and  dynamic  peak  responses  computed  by  the  finite  element 
model  for  the  sphere  are  listed  in  Table  10  (see  Figure  12  for  mesh  loca- 
tions). The  diameter  change  for  the  nonlinear  dynamic  run  is  shown  in 
Figure  17.  Because  the  dynamic  peak  response  is  similar  to  the  static 
peak  response,  only  the  static  results  will  be  discussed. 

The  employment  of  the  nonlinear  models  for  the  tuff  and  the  fiber  rein- 
forced concrete  did  affect  the  structural  response.  All  of  the  relatively 
weak  Tuff  No.  1 elements  became  plastic,  while  all  of  the  concrete  elements 
stayed  elastic.  Thus,  the  difference  between  the  Inear  and  the  nonlinear 
runs  is  entirely  caused  by  the  yielding  of  the  tuff  elti.ients.  The  interior 
vertical  diameter  change  for  the  linear  (-0.32  in)  and  the  nonlinear  (-0.24  in) 
runs  both  exceeded  the  theoretical  change  (-0.17  in).  The  respective  interior 
horizontal  changes  (+0.02  eind  -0.02  in)  were  negligible.  Thus,  the  displaced 
shape  of  the  sphere  at  peak  static  load  is  oval.  This  lack  of  hydrostatic 
behavior  is  also  evident  in  the  radial  stresses  of  the  tuff  elements  adjacent 
to  the  sphere.  For  the  elastic  run,  the  peak  radial  stress  is  -10,500 
psi  at  the  crown  and  -3300  psi  at  the  springline.  For  both  runs,  he  greatest 
concrete  stresses  occur  at  the  springline  interior  in  the  tangential  direction. 
The  Drucker-Prager  tuff  reduces  this  tangential  stress  from  -31,25”  psi 
to  -24,210  psi.  However,  this  is  still  greater  than  the  theoretical  value 
(-15,450  psi).  This  indicates  a negative  arching  behavior  for  the  structure. 
This  behavior  is  also  evident  by  observing  the  vertical  tuff  stresses  at 
the  crown  and  the  springline.  The  stress  at  the  crown  (-9150  psi)  is  greater 
than  the  incident  peak  (-7250  psi),  while  the  stress  at  the  springline 
(-5400  psi)  is  less.  The  non-hydrostatic  loading  produces  both  thrust 
and  moment  in  the  sphere.  Positive  bending  (exterior  fibers  in  compression) 
occurs  at  the  crown  and  invert,  while  negative  bending  occurs  at  the  spring- 
line, The  amount  of  positive  bending  at  the  crown  is  enough  to  cause  the 
highest  compressive  tangential  stresses  to  occur  at  the  exterior  rather 
than  the  interior  surface.  For  the  nonlinear  run,  the  tangential  stress 
is  -13,200  psi  at  the  exterior  and  -7350  psi  at  the  interior. 

Composite  Built-Up  Liner  (Figure  1) 

The  static  and  dynamic  peak  responses  are  listed  in  Table  11.  The 
diameter  change  for  the  nonlinear  dynamic  run  is  shown  in  Figure  18.  Like 
the  homogeneous  sphere,  the  dynamic  peak  response  is  similar  to  the  static 
peak  response.  For  brevity,  only  the  static  results  will  be  discussed. 

Positive  arching  is  indicated  because  the  vertical  tuff  stress  at  the 
crown  is  less  than  at  the  springline.  For  the  nonlinear  run,  the  stress 
at  the  crown  (-3200  psi)  is  less  than  the  incident  peal.  (-7250  psi),  while 
the  stress  at  the  springline  (-9360  psi)  is  greater.  Thus,  the  Elastic 
and  the  Curve  Description  models  which  were  derived  for  the  cellular  concrete, 
act  as  a well-designed  back-packing  material.  The  displaced  shape  for 
both  runs  is  oval.  The  inside  diameter  changes  in  the  vertical  direction 
are  -0.65  in  and  -0.86  in,  and  in  the  horizontal  direction  are  +0.65  in 
and  +.82  in,  for  the  linear  and  nonlinear  runs,  respectively. 
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Unlike  the  Tuff  No.  1 in  the  homogeneous  sphere  mesh,  only  the  tuff 
elements  adjacent  to  the  liner  yield.  This  is  primarily  because  of  the 
increased  strength  of  the  Tuff  No.  2 model.  This  yielding  plus  the  yielding 
in  the  steel  liner  are  the  principal  reasons  for  the  difference  between 
the  linear  and  nonlinear  structural  response.  At  the  maximum  load,  yielding 
of  the  outer  fibers  in  the  steel  elements  near  the  springline  occurred. 
Yielding  is  sufficient  to  cast  doubt  on  the  validity  of  the  linear  run. 
Positive  bending  occurs  at  the  crown  while  negative  bending  occurs  at  the 
springline  of  the  steel  liner.  For  the  nonlinear  run,  the  hoop  stress 
at  the  crown  varies  from  +6800  psi  at  the  inner  surface  to  -44,120  psi 
at  the  outer  surface.  At  the  springline  it  varies  from  -46,220  psi  at 
the  inner  surface  to  -11,430  at  the  outer  surface. 

Composite  Integral  Liner  (Figure  2) 

The  static  and  dynamic  peak  responses  are  listed  in  Table  12.  The 
diameter  change  for  the  nonlinear  dynamic  run  is  shown  in  Figure  19.  Because 
the  dynamic  peak  response  is  similar  to  the  static  peak  response,  only 
the  static  results  will  be  discussed.  The  displaced  shape  of  the  liner 
is  oval.  The  inside  diameter  changes  for  the  linear  and  nonlinear  runs 
in  the  vertical  direction  are  -0.62  in  and  -0.84  in,  and  in  the  horizontal 
direction  are  +0.26  in  and  +0.28  in,  respectively.  The  vertical  tuff  stresses 
at  the  crown  are  -6410  psi  (linear)  and  -5560  psi  (nonlinear),  repsectively. 
This  indicates  positive  arching  behavior.  The  high  strength  of  the  Tuff 
No.  2 model  combined  with  the  structure-media  interaction  produces  no  yielding 
of  the  tuff  elements  in  the  nonlinear  run.  Thus,  the  difference  between 
the  linear  and  nonlinear  results  is  caused  only  by  the  nonlinear  behavior 
of  the  steel  and  concrete.  For  the  nonlinear  run,  the  steel  hoop  stress 
at  the  crown  varies  from  -18,170  psi  at  the  inner  surface  to  -34,420  psi 
at  the  outer  surface.  At  the  springline  it  varies  from  -46,660  psi  to 
-47,860  psi.  Compression  yielding  of  the  steel  liner  occurs  throughout 
the  entire  thickness  from  i4S®  about  the  springline.  This  extensive  yielding 
casts  doubt  on  the  validity  of  the  linear  run  as  is  verified  by  the  pre- 
posterous steel  centroidal  hoop  stress  of  -214,410  psi  at  the  springline. 
Examination  of  the  concrete  hoop  and  radial  stresses  indicates  a thick-wall 
cylinder  behavior  for  the  concrete  shell.  The  radial  stress  at  any  location 
varies  throughout  the  thickness,  while  the  much  larger  hoop  stresses  transmit 
the  load  around  the  cavity.  At  the  maximum  load,  the  concrete  is  plastic 
to  ^45°  from  the  springline  throughout  its  entire  thickness.  The  concrete 
hoop  stress  at  the  crown  varies  from  -5510  psi  at  the  inner  surface  to  -17,840 
psi  at  the  outer  surface,  and  at  the  springline  it  varies  from  -45,570 
psi  to  -28,870  psi.  This  maximum  hoop  stress  at  the  springline  is  approxi- 
mately eight  times  the  unconfined  compressive  strength  (f^  = 5500  psi). 
References  presented  in  Appendix  A indicate  that  plain  concrete  subjected 
to  multi-axial  loading  exhibits  more  ductility  and  higher  strength  than 
concrete  subjected  to  uniaxial  or  biaxial  loading.  The  extremely  high 
concrete  stresses  seem  in  line  with  experimental  data. 
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SU^MARY 


Based  on  these  response  calculations,  the  three  Mighty  Epic  structures 
seem  adequate  for  the  0.5  Kbar  load  level.  However,  this  conclusion  is 
predicated  on  the  appropriateness  of  the  material  models.  This  is  of  par- 
ticular concern  for  the  concrete  models,  which  were  derimed  from  little 
more  than  estimated  values  of  f^.  Significant  confidence  in  calculations 
can  not  be  established  until  better  material  data  is  gathered  and  its  validity 
checked  against  prototype  structures. 
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'Location  of  response,  see  Figure  12.  All  valuer  nf  siress/stratn  taken  from  elemeni  centroid. 

^Node  number  asaoctated  sirirh  dtspUceirtMit  {inches)  (6^  verticai.  and  iy,  horizontal) 

'Element  number  associated  srith  streas/strain  ia  and  f vertical,  o and  e horizontal,  and  a arul  ihear)  (psi/M  in  per  in). 


) III  i-'piO  * 


'tocitkHi  of  reiponK.  tec  FlRure  12.  All  valuet  of  atrctt/itriin  Uken  from  element  centroid. 

’Node  tttimber  eaDdeted  vKh  doptecemenii  (inches)  (6^  mtlctl,  end  iu  borizontti). 

'EJetnent  number  taociated  wHh  strea/ttrsin  (a  . and  c mtkal.  Oy  and  ry  horizontal,  arul  0 y aid  «,y  daar)  (ptl/>i  bi  per  in). 
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■Location  of  retpontc.  Fifurr  1 2 All  valuer  of  ttma/airain  taken  from  element  centroid. 

^Node  number  aaaociated  with  diaplaccmenK  (Ay.  vertical,  and  A||.  horiiontal). 

<FJemetH  number  ataocMtcd  with  •treWetrain  (o^  and  wrtkal,  oy  •'*<1  <Y*  f>onzontai.  and  o^y  and  c^y,  drear). 

<k>nly  the  peak  reaponiea  arc  drown.  The  peak  dHpiacementt  occur  at  the  end  of  the  run  at  17.6  tnaec  while  the  peak  atreaBra  and  atraiira  occur  at  the  timea  iirdkatrd  in  parenthevea. 
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COMPOSITE  BUILT-UP  LINING 
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Figure  1.  Homogeneous  Sphere  and  Composite  Built-up  Liner 


(a)  Free  Field 


(b)  Other  Structures 

Dj  is  the  of  the  region  of 

interest  - 

Dg  Is  the  additional  depth  of  mesh 
necessary  to  avoid  reflections 


Figure  3.  General  Character  of  Finite  Element  Meshes 


Dovmstream 


Figure  4.  Predicted  Free  Field  Radial  Stress  Time  Histories  for  Mighty  Epic  (Uniaxial 


Figure  5.  Predicted  Free  Field  Velocity  Time  Histories  for 
Mighty  Epic  (uniaxial  case.  Ref.  3) 
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Figure  8.  Free  Field  Response  Stress  vs.  Time  (Element  24) 
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fb)  CASE  1,  control 
c^se  made  of  all 
4 node  elements 


(c)  CASE  2,  transition 
of  4 node  to  8 node 
elements 


(d)  CASE  3,  transition 
of  8 node  to  4 node 
elements 
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Figure  10.  Meshes  used  to  check  element  compatibility 


(a)  Mesh  3.0 


Figure  11.  Meshes  for  Free  Field  with  Cylindrical  Cavity  (Note:  thin  cylindrical 

cavity  liner  too  small 

to  be  seen  on  figures) 
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1.  All  stresses  are  at  the  centroids  of  the  elements  nearest  the  locations  shovm. 

2.  A-I  are  either  in  tuff  or  at  tuff /liner  interface  where  A is  top  left.  B top  mid, 
C is  top  right,  D at  crown,  E adjacent  springline,  F adjacent  E,  G at  mid  mesh, 

H at  right  edge  even  with  springline,  I at  invert. 

3.  J,  K,  and  L are  at  the  interior  surface  of  the  liner  where  J is  at  crown,  K at 
springline,  and  L at  invert. 


Figure  12.  Locations  of  Responses  shown  in  Tables  2-12 
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Change  in  Diameter  (in) 


Horizontal 


Figure  13.  Change  in  diameter  history-Run  DCS 
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a.  Tuff  Elements  (0"  ~ 504") 

Figure  16.  Composite  Integral  Liner  (Mesh  6.0) 
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APPENDIX  A 

DEVELOPMENT  OF  MATERIAL  PARAMETERS 


INTRODUCTION 

The  material  models  employed  should  accurately  represent  the  actual 
materials  under  field  conditions.  That  is,  the  constitutive  relationships 
must  be  valid  throughout  the  expected  loading  range.  Several  linear  and 
nonlinear  models  are  available  in  NONSAP  to  describe  material  behavior: 
the  linear  elastic  isotropic  model,  the  von  Mises  model,  the  Drucker-Prager 
model  and  a curve  description  model  with  tension  weakening.  The  von  Mises 
and  Drucker-Prager  models  have  explicit  yield  functions  which  specify  the 
state  of  multiaxial  stress  corresponding  to  the  start  of  plastic  flow. 

The  von  Mises  yield  surface  is  parallel  to  the  hydrostat,  and  the  Drucker- 
Prager  yield  surface  is  at  an  angle  to  it.  Both  models  may  incorporate 
flow  rules  and  hardening  rules.  Currently  in  NONSAP  the  von  Mises  condition 
can  be  used  for  either  perfectly  plastic  conditions  or  isotropic  hardening; 
and  the  Drucker-Prager  yield  condition  is  limited  to  elastic,  perfectly 
plastic  analysis.  In  both  cases  an  initially  elastic  material  is  assumed. 

The  curve  description  model  is  an  incremental  stress-strain  law  which  de- 
scribes tangent  bulk  moduli  in  loading  and  unloading  (Ky)  and  the 

loading  shear  modulus  (Gl)  as  piecewise  linear  functions  of  volumetric 
strain  (Oy) . TTiis  model  has  no  explicit  yield  function,  and  it  allows 
the  material  to  weaken  abruptly  (crack)  if  the  tensile  stress  resulting 
from  applied  loads  exceeds  the  allowable  tensile  strength.  Model  selection 
for  the  analysis  was  based  on  laboratory  material  test  data  supplied  by 
Defense  Nuclear  Agency  (DNA)  contractors  and  supplemented  with  material 
behavior  information  from  the  technical  literature. 

TUFF 

A1 

The  experimental  test  data  supplied  by  Paul  Weidlinger  Associates 
was  initially  used  to  derive  a tuff  model.  Results  from  one-dimensional, 
hydrostatic,  and  triaxial  tests  are  shown  in  Figures  Al,  A2  and  A3,  respec- 
tively. These  figures  indicate  that  the  material  does  not  behave  linear 
elastically  at  the  stress  levels  of  interest  and  that  some  form  of  yielding 
must  be  incorporated  into  the  constitutive  relationships.  The  Drucker-Prager 
yield  model  was  selected  to  model  this  inelastic  behavior.  The  free  field, 
the  cavity,  and  the  homogeneous  sphere  used  this  model  (Tuff  No.  1).  The 
composite  built-up  liner  and  the  composite  integral  liner  used  a different 
tuff  model  (Tuff  No.  2)  based  on  more  recent  experimental  data  from  Terra 
Tek^.  The  derivation  of  Tuff  No.  1 is  outlined  below.  A similar  techniqae, 
which  is  not  shown,  was  also  used  for  Tuff  No.  2. 

Mohr  circles  are  shown  in  Figure  A4  for  the  unconfined  compression 
test  and  the  0.5-Kbar  confining  stress  test.  A failure  envelope  was  then 
constructed  which  has  a cohesive  intercept,  c = 0.3  Kbar  and  a friction 
angle,  ((>  = 9.23°.  Values  of  E and  v were  chosen  to  best  describe  the  behavior 
of  the  tuff  at  the  0.5-Kbar  stress  level.  A trial  and  error  technique 
is  used  to  determine  these  values  from  the  given  data  and  from  the  following 
elastic  equations: 
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where: 


Cl  = Modulus  of  Elasticity  in  shear  (psi) 

K = Bulk  Modulus  (psi) 

M = One-dimensional  confined  modulus  (psi) 

K^=  At-rest  coefficient  of  lateral  earth  pressure 

Because  1-D  strain  tests  are  more  difficult  to  perform  on  rock  than 
soil,  the  data  from  the  hydrostatic  and  triaxial  tests  are  used  to  determine 
E and  v.  The  1-D  test  data  are  used  to  check  results.  Values  for  the 
shear  modulus,  G,  at  each  level  of  confining  stress  are  calculated  from 
Figure  A3  (G  = one-half  the  slope  of  stress  difference  versus  strain  dif- 
ference) . These  values  were  then  plotted  against  the  confining  stress 
in  Figure  AS.  The  trial  and  error  procedure  is  as  follows: 

1.  Assume  a value  for  K^. 

2.  Determine  the  hydrostatic  stress  state  equivalent  to  the  0.5  Kbar 
1-D  confined  stress  state. 


where:  a = Vertical  1-D  stress 

V 

Pj^l  = Equivalent  hydrostatic  stress 

3.  Determine  the  secant  bulk  modulus,  K,  at  this  stress  level  from 
Figure  A2. 

4.  Determine  the  confining  stress,  pj,  in  the  triaxial  test  equivalent 
to  the  0.5-Kbar  1-D  confined  stress  state. 

Pt  = *^0  % 

5.  Determine  G at  this  stress  level  from  Figure  AS. 

6.  Determine  E and  v by  solving  Equations  1 and  2 simultaneously. 
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7.  Determine  Kq  from  Equation  4 and  compare  it  to  the  assumed  value. 


Iterate  these  7 steps  until  the  results  converge.  Assuming  an  initial 
value  of  Kq  = 0.30,  the  following  results  are  obtained  after  three  iter- 
ations : 

Kq=  0.40 

G = 12  Kbar  = 174,000  psi 

K = 24  Kbar  = 348,000  psi 

E = 30.9  Kbar  = 448,000  psi 

V = 0.286 

M = 40  Kbar  = 581,000  psi 

The  dotted  straight  lines  in  Figures  Al,  A2  and  A3  are  drawn  to  repre- 
sent these  moduli.  The  value  of  M obtained  directly  from  Figure  Al  at 
the  0.5-Kbar  stress  level  is  32  Kbar.  This  compares  reasonably  well  with 
the  value  calculated  above. 

These  derived  values  for  Tuff  No.  1,  as  well  as  the  Tuff  No.  2 values, 
are  shown  in  Table  1, 

STEEL 

The  cavity,  the  composite  built-up  liner,  and  the  composite  integral 
liner  contain  a steel  cylindrical  liner.  The  cavity  liner  is  assumed  to 
behave  elastically  while  the  von  Mises  model  is  selected  to  describe  the 
nonlinear  behavior  of  the  other  two.  The  material  parameters  are  shown 
in  Table  1. 

FIBER  REINFORCED  CONCRETE 

The  homogeneous  spherical  shell  is  constructed  of  steel  fiber  reinforced 
concrete.  In  addition  to  increased  tensile  strength,  concrete  containing 
steel  fibers  has  slightly  increased  compressive  strengths,  greater  strain 
capacity,  toughness,  post-peak  strength,  and  post-peak  integrity  than  con- 
ventional concretes. These  characteristics  are  further  enhanced  under 
biaxial  and  triaxial  compression.^^  Thus,  the  material  definitely  possesses 
elasto-plastic  characteristics. 

The  designers  of  the  structure  have  indicated  that  a shot  day  strength 
of  17,000  psi  is  probable.  However,  concrete  experts  at  CEL  feel  that 
the  interaction  and  compounding  of  the  various  factors  contributing  to 
such  a high  strength  are  not  well  understood.  A more  reasonable  assumption 
for  the  shot  day  dynamic  unconfined  compressive  strength  is  10,000  psi*. 

This  value  is  used  in  the  analysis. 

* six  months  a^ter  shot  day,  CEL  tested  concrete  cores  from  an  unused  sphere. 
The  average  static  compressive  strength  was  12,450  psi. 
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A technical  literature  search  dealing  with  triaxial  testing  of  concrete 
was  carried  out  to  select  the  appropriate  plasticity  model  for  the  fiber 
reinforced  concrete.  Literature  indicates  that  concrete  behavior  is  similar 
to  tuff.  That  is,  the  Mohr  failure  envelope  has  both  a cohesion  intercept 
and  a friction  angle.  Thus,  the  Drucker-Prager  model  is  used.  Since  triaxial 
test  data  for  the  fiber  reinforced  concrete  was  not  available,  c and  (ji 
were  estimated  from  the  data  available  in  the  literature.  Triaxial  test 
data  from  References  AS  and  A6  are  used  to  derive  values  of  c and  That 
is,  Mohr*s  circles  were  plotted  using  principal  stress  differences  normalized 
with  respect  to  f^  and  the  average  failure  envelope  is  drawn.  The  results 
indicate  that  the  cohesion  is  approximately  0.36  f^  and  the  friction  angle 
approximately  35°.  Thus,  for  10,000  psi  fiber  reinforced  concrete,  the 
cohesion  is  3600  psi.  According  to  Reference  A3,  the  addition  of  steel 
fibers  to  concrete  does  not  significantly  affect  Young's  modulus  or  Poisson's 
ratio.  An  average  value  for  Poisson's  ratio  of  0.24  was  selected  from 
Reference  A3.  The  following  equation  was  used  to  determine  Young's  modulus: 

E = 33  w^’^'Vf^  (5) 

where,  w = weight  density  of  concrete  (Ib/ft^) . For  144  Ib/ft^  concrete, 
the  calculated  E is  5,700,000  psi. 

REINFORCED  CONCRETE 

The  cylindrical  composite  integral  liner  consists  of  a thin  steel 
liner  surrounded  by  reinforced  concrete.  The  presence  of  the  various 
rebars  in  the  plain  concrete  shell  changes  the  outer  shell  from  a relatively 
simple  homogeneous  structure  to  a complicated  "matrix"  of  interacting 
materials.  The  effect  of  the  rebar  on  the  structural  behavior  must  therefore 
be  considered  in  the  analysis. 

A special  finite  element  study  of  a concentrated  load  acting  on  the 
end  of  a cantilever  beam  was  undertaken  to  evaluate  the  interaction  between 
the  concrete  and  rebar.  TTie  three  meshes  shown  in  Figure  A6  were  used 
to  model  the  beam,  which  was  similar  in  thickness  (5  in)  to  the  reinforced 
concrete  shell.  Meshes  1 and  2 represent  a plain  concrete  beam;  Mesh  3 
represents  a reinforced  beam.  The  concrete  elements  are  all  8-node  quad- 
rilaterals; the  steel  rebar  is  modeled  by  truss  elements.  The  area  of 
the  truss  members  was  adjusted  to  match  the  steel  area  of  the  actual  rein- 
forced concrete  shell.  The  primary  objective  is  to  determine  the  influence 
of  the  rebar  on  structural  response.  A concentrated  load  (P  = 1,000  lb) 
was  applied  to  the  beam  tip.  The  theoretical  end  deflections  of  the  plain 
concrete  and  the  reinforced  concrete  beams  were  0.117  inch.  The  deflections 
computed  by  NONSAP  were  0.120  inch  for  all  three  meshes.  This  indicates 
that  a plain  concrete  mesh  can  be  used  to  model  the  stiffness  of  a reinforced 
concrete  shell. 

Under  the  postulated  loading  conditions,  a triaxial  compressive  stress 
state  should  exist  throughout  the  reinforced  concrete.  Recent  experiments^S, 
A6,  A7  indicate  that  plain  concrete  subjected  to  multiaxial  loading,  exhibits 
more  ductility  and  higher  strength  than  concrete  subjected  to  uniaxial 
or  biaxial  loading.  This  behavior  suggests  that  a Drucker-Prager  plasticity 
model  is  applicable. 
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The  only  material  property  data  supplied  to  CEL  was  the  unconfined 
compressive  strength  (f^)  of  the  plain  concrete,  5500  psi.  The  previous 
approach  of  determining  the  friction  angle  and  cohesion  for  the  fiber  rein- 
forced concrete  is  also  used  for  the  plain  concrete.  The  friction  angle 
equals  35°  and  the  cohesion  is  0.36  f'  (1980  psi).  The  value  of  E calculated 
from  Equation  5 is  4,264,000  psi,  while  the  value  of  Poisson's  ratio  is 

0.17. 

CELLULAR  CONCRETE 

Cellular  concrete  is  a lightweight,  low  modxlus  backpacking  material 
surrounding  the  steel  liner  in  the  composite  built-up  liner.  Its  function 
is  to  reduce  the  stresses  acting  on  the  structure  through  positime  arching 
action. 

Laboratory  test  data  on  cellular  concrete  performed  by  Waterways  Experi- 
mental Station  was  used  to  determine  which  NONSAP  material  model  is  appli- 
cable. The  idealized  results  from  a hydrostatic  test  on  60  pcf  cellular 
concrete  are  shown  in  Figure  A7.  Three  distinct  regions  of  behavior  are 
present.  Initially,  the  material  deforms  linearly  up  to  a yield  stress 
(Oy^)  of  about  1500  psi.  At  this  stress  level  the  material  begins  to  crush, 
reducing  the  amount  of  initial  air  space.  This  ideal  plastic  behavior 
stops  at  a hardening  strain  (Cj^.)  of  about  20  to  30  percent.  At  this  strain 
value,  the  void  space  has  essentially  been  reduced  to  zero.  The  material 
is  now  capable  of  resisting  subsequent  loading.  The  modulus  in  this  region 
is  approximately  one-half  the  initial  modulus.  This  material  behavior 
can  best  be  approximated  in  NONSAP  by  the  Curve  Description  Model.  Neither 
the  Drucker-Prager  nor  the  von  Mises  model  is  capable  of  simulating  this 
unique  elastic-plastic  behavior.  To  obtain  the  input  for  the  Curve  Descrip- 
tion Model,  the  following  trial  and  error  procedure  is  used: 

1.  Determine  approximate  values  for  the  model  based  on  Figure  A7. 

2.  Hydrostatically  load  an  axisymmetric  element  using  NONSAP. 

3.  Plot  a p-ey  curve  and  compare  with  Figure  A7. 

4.  Modify  the  model  and  repeat  steps  2 and  3 until  satisfactory  con- 
vergence is  obtained. 

The  final  values  for  the  model  are  outlined  in  Table  Al.  The  NONSAP  p-e^ 
curve  for  this  model  is  shown  in  Figure  A7. 

The  Elastic  Model  shown  in  Table  1 was  obtained  by  assuming  a Poisson's 
Ratio  of  0.48  and  a Bulk  Modulus  of  7500  psi  (see  Figure  A7) . Like  the 
actual  cellular  concrete  behavior,  this  linear  approximation  is  very  soft 
and  has  little  shear  resistance. 
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Figure  A1 . Results  of  static  one-dimensional  compression  test 
on  tuff 
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Mesh  3 


Concrete  Elements 
□ Steel  Elements 


Figure  A6.  Cantilever  beam  meshes. 


Figure  A7.  Hydrostatic  test  results 


APPENDIX  B 

ELEMENT  AND  STEP  SIZE  SELECTION 


In  order  to  theoretically  determine  the  element  and  step  size  of  the 
columnar  mesh  (shown  in  Figure  Bla) , two  analogous  systems  are  employed 
(Figure  Bib  and  Blc) . Bib  is  a lumped  mass  and  spring  system  defined  as 
follows : 

k = ^ Constrained  spring  stiffness  (1) 

1j 

m = pAL  Lumped  Mass  (2) 


where : 


A is  the  element's  cross  sectional  area 
M is  the  tuff's  constrained  modulus 
L is  the  element's  length 
p is  the  tuff's  mass  density 

Blc  is  an  ideal  continuous  constrained  rod  system  with  the  following  wave 
speed  for  the  tuff: 

^ ”s]^=  57.5  in/msec  (3) 

To  enable  the  spring/maso  system  to  reasonably  approximate  the  rod's  response 
to  a ramp  loading,  P(t),  the  following  expression  derived  from  Reference 
B1  may  be  employed: 


where : 

t^  is  the  rise  time  of  P(t)  (defined  in  Figure  6b) 

T is  the  period  of  a single  spring/mass 

Inverting  this  expression,  the  maximum  frequency  (fjnax)  which  needs  to 
be  transmitted  through  the  spring/mass  system,  given  a rise  time  tj.,  is: 

f = 3 = 3 = 500  cps  (4) 

t 006 
r 

The  maximum  frequency  that  can  be  propagated  through  the  spring/mass  system 
is  given  in  Reference  B2  as  follows: 
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_c_ 

Ltt 


(5) 


Substituting  from  (1) , (2) , and  (3) : 


f = 


_1_ 

Ltt 


TTie  maximum  frequency  Cfninv)  that  can  be  expected  to  propagate  unattentuated 
through  system  (Bib)  is  given  by  (5).  Combining  (4)  and  (5)  yields: 


3 


= f _ 


max 


_C_ 


(6) 


L 


= iil  ^ (57.5)  (6)  ^ 
3Tr  3n 


36.6  inches 


An  alternative  mesh  size  selection  procedure  is  given  in  Reference 
B3  for  use  in  conjunction  with  calculations  for  the  continuous  rod  system 
as  follows: 


L 


I 


L = 6(57.5)  =115  inches 


To  select  a time  step.  At,  for  use  in  conjunction  with  the  implicit 
step-by-step  integration  procedure  of  the  NONSAP  program  the  following 
criteria  are  derived  for  a single-degree-of-freedom  system  (Reference  B4) : 


Components  of  P(t)  with  periods  less  that  5At  will  experience 

attenuation  caused  by  the  time  integration  algorithm.  Thus,  to 

pass  f requires : 

^ max 


At  < 


1 


?7 


max 


1 = 0.4  msec 

5T5(5ST 


The  wave  propagation  through  three  columnar  meshes  (Figure  B2)  was 
studied  to  determine  the  appropriateness  of  the  element  and  step  size 
equations.  The  dynamic  loading  shown  in  Figure  6b  was  applied  to  the  top 
surface  of  each  mesh  (linear  Tuff  No.  1).  The  vertical  stress  histories 
at  the  locations  shown  in  Figure  B2  are  plotted  in  Figure  B3.  Figure  B4 
shows  the  shape  of  the  stress  wave  as  a function  if  L and  At.  Tbe  data 
shown  indicates  that  for  L = 40  inches  and  At  = 0.4  msec  satisfactory  results 
are  achieved. 
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(a) 

lamped  mass 
Finite  element 
system 


(b)  (c) 

lumped/mass-  continuous 

constrained  constrained  rod 

spring  system 


where 

M = 581,380  ps 
P = .00017577 


{constrained  modulus  for 
the  tuff,  see  Table  1. 
lb  -560“  mass  density  of  tuff. 


0.286 


Poisson's  ratio  for  tuff 


Figure  Bl.  Columnar  Systems  Used  For  Tuff  Free  Field  Calculations 
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Time  (msec) 

(a)  190-inch  depth  (Element  10)  , L=20  inch 
Figure  B3.  Vertical  Stress  vs.  time  for  Free  Field  Meshes 


40  inches,  At  = 0.2,  0.4,  and  1.0msec 


APPENDIX  C 

STUDY  OF  COMPATIBILITY  OF  4-\'ODE  AND 
HIGHER  ORDER  QUADRILATERAL  ELEMENTS 


In  some  circumstances  it  is  advantageous  to  join  8-node  to  4-node 
quadrilateral  elements^  Three  different  column  meshes  are  employed  to 
investigate  the  effect  on  stresses  and  displacements  (Figure  10).  Case  1 
represents  the  "ideal"  mesh,  containing  two  columns  of  4-node  elements. 

The  Case  2 mesh  consists  of  two  columns  of  4-node  elements  in  the  upper 
200  inches  while  the  lower  600  inches  contain  one  column  of  8-node  elements. 
The  Case  3 mesh  consists  of  one  column  of  8-node  elements  in  the  upper 
200  inches  with  the  lower  600  inches  consisting  of  two  columns  of  4-node 
elements.  For  all  elements,  L is  20  inches.  A O.SOKbar  ramp  loading 
(t_  = 6 msec)  is  applied  to  the  top  surface  of  each  mesh.  Both  linear 
and  Drucker-Prager  Tuff  No.  1 are  used  for  the  material  characterization. 

The  response  at  the  200-inch  depth  interface  and  just  above  and  below 
is  investigated.  Each  run  is-for  30  time  steps  where  At  = 0.4  msec.  The 
reported  stresses  and  displacements  occur  at  t = 12  msec  (30th  time  step); 
these  responses  represent  the  maximums  expected  for  this  region.  Tables  Cl 
and  C2  show  the  responses  at  locations  indicated  by  Figure  Cl  using  an 
integration  order  of  three.  To  measure  the  impact  of  integration  order 
on  the  response,  the  stresses  at  all  integration  points  are  shown  in  Tables  C3 
and  C4  for  integration  orders  2 and  3. 

The  results  indicated  that  for  integration  order  3,  such  mesh  transitions 
produce : 

1.  No  significant  change  in  overall  displacements. 

2.  Element  stresses  that  are  not  significantly  affected  beyond  two 
elements  of  the  transition. 

3.  Correct  centroidal  stresses  for  only  the  4-node  quad  elements 
near  the  region.  Also,  the  average  stress  of  all  nine  stress  points  is 
correct  for  only  the  4-node  quad  elements. 

4.  Neither  centroidal  stress  nor  the  average  stress  for  the  8-node 
quad  elements  are  correct  near  the  transition. 

For  integration  order  2: 

1.  No  significant  change  in  overall  displacements. 

2.  Element  stresses  that  are  not  significantly  affected  beyond  two 
elements  of  the  transition. 

3.  Near  the  interface,  all  four  stresses  for  the  8-node  quad  elements 
are  correct. 

4.  Near  the  interface,  all  four  stresses  for  the  4-node  quad  elements 
are  incorrect.  However,  the  average  stress  is  correct. 


Table  C3.  Vertical  Stress  and  its  Deviation^  Caused  by  Transition  for 
Linear  Tuff  (Integration  Order  2)  (psi) 


Stress*’ 

Points 

Location  of  Stress^  for 

Element  No.  (Case  2) 

9 

19 

10 

20 

21 

22 

1 

-7308 

-61 

-7207 

40 

-6784 

407 

-7611 

-420 

-7190 

-10 

B 

2 

-7308 

-61 

-6784 

407 

1 

-7611  ' 

-420 

1 

B 

-7249 

-1 

3 

-7308 

-61 

-7611 

-420 

-6784 

407 

-7190 

-10 

B 

4 

-7207 

40 

m 

-7611 

-420 

-6784 

407 

-7184 

-A 

-7249 

-1 

Element  No.  (Case  3) 

9 

10 

11 

41 

12 

42 

1 

-7225 

22 

-7177 

14 

-6780 

400 

-7607 

-All 

-7332 

-84 

m 

i 

2 

m 

-7181 

10 

-6780 

400 

-7607 

-427 

-7332 

-84 

B 

3 

-7225 

22 

-i\n 

14 

-6780 

400 

B 

B 

4 

1 

B 

-7181 

10 

-7607 

^27 

-7232 

16 

m 

‘‘First  entry  in  box  is  vertical  stress  while  the  second  entry  is  the  difference  between 
first  entry  and  corresponding  stress  for  the  “ideal"  solution  (Case  1,  Figure  10b) 
^Location  of  stress  output,  4 per  element,  see  Figure  10a. 

‘"Element  numbers  shown  in  Figure  10c  and  lOd. 
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ruble  C'4.  Verlieul  Sliess  uikI  i(s  IX’vialion^  Caused  by  Transition  (or 
Linear  Tult  (Integration  Order  3)  (psi) 


Loeation  of  Stress^  for 


■) 

■D 

-7519 

-6704 

-7218 

mm 

-32S 

476 

■■ 

■1 

30 

1 

•} 

-7 1 ‘)4 

-7522 

-6704 

-7677 

-7218 

53 

-33 1 

476 

-4<n 

30 

continued 
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Table  C4.  Continued 


f 


Stress*’ 

Point 

Location  of  Stress'"  for 

Element  No.  (Case  3) 

9 

10 

1 1 

41 

12 

42 

4 

-7287 

-40 

-6756 

435 

-7190 

-10 

-7190 

-10 

-7277 

-29 

B 

5 

-7296 

-49 

-7190 

-10 

-7190 

-10 

-7277 

-29 

-7277 

-29 

b 

m 

m 

-7190 

-10 

-7190 

-.0 

-7277 

-29 

-7277 

-29 

7 

1 

-7176 

71 

-7515 

-324 

-6704 

476 

-7218 

30 

B 

8 

m 

-7519 

-328 

-7677 

-497 

-6704 

476 

-7218 

30 

B 

0 

-7194 

53 

-7522 

-331 

-7677 

^97 

W 

-7218 

30 

B 

“First  entry  in  box  is  vertical  stress  while  the  second  entry  is  the  difference  between 
first  entry  and  corresponding  stress  for  the  “ideal”  solution  (Case  1,  Figure  lOb) 
^Location  of  stress  output,  9 per  element,  see  Figure  10a. 

‘■Element  numbers  shown  in  Figure  10c  and  lOd. 
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APPENDIX  D 

RESTART  USING  PROGRESSIVELY  LARGER  At 


The  columnar  mesh  shown  in  Figure  D1  is  used  to  demonstrate  that  pre- 
dictions at  the  structure  site  can  be  satisfactorily  computed  using  a pro- 
gressively increasing  At.  The  Mighty  Epic  forcing  function  (Figure  6a) 
is  applied  to  the  top  surface  of  the  linear  Tuff  No.  1 column.  For  this 
run: 

1.  At  = 0.4  msec  for  the  first  50  steps  (20  msec  total) 

2.  At  = 1.0  msec  for  next  50  steps  (70  msec  total) 

3.  At  = 5.0  msec  for  next  six  steps  (100  msec  total) 

4.  At  = 50  msec  for  the  next  28  steps  (1500  msec  total) 

The  stress  history  at  the  290-inch  depth  (Element  15)  is  shown  in  Figure  D2. 

For  this  idealization  the  reflected  wave  arrives  at  Element  15  in  204  msec. 

As  shown  in  Figure  D2,  the  amplitude  of  the  reflected  stress  wave  is  reduced 
by  9S%. 

This  study  suggests  that  a progressively  increasing  At  is  a viable 
means  of  reducing  the  run's  cost  and  the  system's  degrees-of-freedcxn  for 
"long-time"  finite  element  runs.  However,  this  conclusion  must  be  checked 
on  a "case-by-case"  basis  with  verification  runs  being  made  for  the  particular 
idealization  involved.  Special  attention  must  be  applied  in  the  region 
of  interest  to  insure  that  secondary  transient  responses  (i.e.,  those  created 
by  the  structure's  response  to  the  forcing  function)  are  neither  prematurely 
attenuated  nor  wastefully  computed. 
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APPENDIX  E 

USE  OF  VELOCITY  AS  A FORCING  FUNCTION 


Tlie  velocity  time  history  at  the  structure  location  ("0"  on  Figure  5) 
is  integrated  yielding  a displacement  time  history.  Stiff  springs  are 
inserted  at  the  top  nodes  of  a columnar  mesh  (Figure  El) . Loads  (Figure  E2) 
which  cause  the  desired  displacement  time  history  are  applied  to  these 
nodes.  Both  linear  and  Drucker-Prager  Models  are  used  to  describe  the 
material.  The  vertical  stress  vs  time  is  shown  in  Figures  E3  and  E4  and 
compared  with  the  given  radial  stress  time  histories  (i.e.,  "0"  wave  of 
Figure  4).  Both  models  show  good  agreement  up  to  the  peak  stress  level. 
Niether  model  agrees  after  this  time.  Note:  When  integration  order  3 
was  used,  indefinite  results  were  detected.  Only  integration  order  2 could 
be  used. 
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Figure  E3.  Radial  (vertical)  stress  in  Element  1;  linear  elastic  Tuff 
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